# necessite les packages : car et gap
# necessite les scripts GWA_analysis_functions.R et qqmultiple.R

# chargement des fonctions
source("GWA_analysis_functions.R")

################# VUE D'ENSEMBLE DES DONNEES et elimination de SNPs #####################

# recuperation de l'ensemble des donnees
create_summary_data(inbfile="GWAdata", del_AT_GC=F, del_heterosome=F, outputfile="summary_data.txt")

#creation d'une liste de SNPs à eliminer suivant Hardy-Weinberg, pourcentage de missing et maf
snp_selection(HW=0.0001, miss=0.05, maf=0.1,summary_file="summary_data.txt")# outputfile = "exclusionSNP.txt"

# elimination des SNPs en fonction des F_U et F_A
del_freq(inbfile="forStrat", th_freq=0.05, outputfile="SNPoutfreq.txt")


########################## STRATIFICATION ############################

# fonction permetttant de faire le plot de tous les PI HAT, le plot des moyennes de PI HAT par individu et l'histogramme de la distribution des moyennes de PIHAT 
PIHAT.plots(genomefile="plink.genome", general_plot_title="PI-HAT", mean_plot_title="average PI-HAT for each individual", distrib_title="distribution", distrib_breaks=100, outputfile="plots_PI-HAT")

# creation d'une liste d'individu à exclure de l'analyse
PIHAT.exclusionlist(genomefile="plink.genome",famfile="forStrat2.fam",th_indiv=0.1,th_mean=0.015,controls=c("ADN103M2","ADN103M3","ADN103M4","ADN103M5","ADN103M6"))

# creation d'un fichier phenotype pour les composantes
create_strat_phenofile(mdsfile="strat.mds", outputfile="pheno_strat.txt")

# creation de plots montrant l'association des individus aux composantes
plot.indiv.components(nb_comp=4, famfile="forStrat2.fam", phenofile="pheno_strat.txt", outputfile="qqplot_components", color="sex")#color=sex for coloring by sex, color=status for coloring cases/controls

# creation de plots montrant l'association d'individus aux composantes en fonction de l'etude desquels ils sont tires
plot.indiv.components.orig <- function(nb_comp=4, all_famfile="all.fam", casesctrls_famfile="pvmcont.fam", ctrls_famfiles=c("oba.fam","canada.fam"),phenofile="pheno_strat.txt",outputfile="qqplot_components_orig")

# creation d'un unique plot avec une ligne horizontale ou verticale
plot.indiv.components.abline(comp_number=c(1,3),famfile="forStrat2.fam", phenofile="pheno_strat.txt",color="status",Hline=0.028,Vline=NULL)#Hline : ligne horizontale, Vline : ligne verticale

# creation d'un tableau bilan des donnees d'association aux composantes
summary.qassoc(nb_comp=4, assocfile="pheno", outputfile="summary-qassoc.txt")

# creation de qqplot d'association aux composantes
qqplot.components(nb_comp=4, IC=0.95, summaryfile="summary-qassoc.txt", outputfile_sq="qqplot_statsq_assoc_components", outputfile_p="qqplot_pval_assoc_components")


######################### ANALYSE PAR PLAQUE #############################
###(creation des qqplots de comparaison de plaque les unes contre les autres et manhattan plot par plaque)###

# creation des fichiers cluster - a faire pour les cas et les controles
create_cluster_file(inbfile="data", plate_file="Plaks.txt", num_Pdel=NULL, status=1) # num_Pdel : numero de plaque à supprimer, status : 1 for controls, 2 for cases

# creation du fichier de phenotype - a faire pour les cas et les controles
create_pheno_file(cluster_file="cluster-ctrls.dat", status=1) #status : 1 for controls, 2 for cases

# creation des fichiers d'association sous plink
association()

#creation d'un tableau d'informations sur les chi² et les p-values
info_chisq_pval(p.names="P", status=1)

# creation d'un QQplot par plaque - a faire pour les cas et pour les controles
qqplot.p(p.names=c("P1","P2","P3"), nb_by_graph=4, IC=0.95, status=2, outputfile="qqplot") #nb_by_graph : nb de plaque par graphique ; status : 1 for controls, 2 for cases

# creation d'un manhattan plot pour chaque plaque
mhtplot.p(p.name="P1", status=2, title="manhattan plot", outputfile="mhtplot1") #status : 1 for controls, 2 for cases


################################## EXAMPLES #####################################

## une stratification
# chargement des fonctions
source("GWA_analysis_functions.R")
# elimination des SNPs en fonction des F_U et F_A
del_freq(inbfile="forStrat", th_freq=0.05, outputfile="SNPoutfreq.txt")
# recuperation de l'ensemble des donnees (ici on supprime ce qui correspond aux chromosomes X, Y et MT ainsi que les couples d'allèles du type : A/T, T/A, G/C, C/G)
create_summary_data(inbfile="GWAdata", del_AT_GC=T, del_heterosome=T, outputfile="summary_data.txt")
# Suppression des SNPs selon leur equilibre de hardy-weinberg, maf, et pourcentage de missing
snp_selection(HW=0.0001,miss=0.05,maf=0.1)
inbfile <- "GWAdata"
out <- "forStrat"
system(paste("plink --noweb --bfile ",inbfile," --exclude exclusionSNP.txt --make-bed --out ",out,sep=""))
# Suppression des SNPs en fonction de leur frequence chez les cas et les contrôles 
del_freq(inbfile="forStrat",th_freq=0.05)
system(paste("plink --noweb --bfile ",inbfile," --exclude SNPoutfreq.txt --make-bed --out ",out,sep=""))
# Selection des SNPs independants
system(paste("plink --noweb --bfile ",inbfile," --hwe 1e-04 --geno 0.05 --maf 0.1 --mind 0.05 --indep-pairwise 50 10 0.2",sep=""))
system(paste("plink --noweb --bfile ",inbfile," --extract plink.prune.in --make-bed --out ",out,sep=""))
# Elimination des SNPs en LD residuel
bfile <- out
system(paste("plink --noweb --bfile ",bfile,"--ld-window 50 --ld-window-kb 5000 --ld-window-r2 0.2 --out ld",sep=""))
threshold <- 0.19999
system(paste("perl snpIn.pl ld.ld",threshold))
out <- "forStrat2"
system(paste("plink --noweb --bfile ",bfile,"--exclude SNP_out.txt --make-bed --out ",out,sep=""))
bfile <- out
# Determination de la matrice IBD
system(paste("nohup plink --noweb --bfile ",bfile,"--genome &",sep="")) 
# plots PIHAT
PIHAT.plots()
PIHAT.exclusionlist(genomefile="plink.genome",famfile="forStrat2.fam",th_indiv=0.1,th_mean=0.015,controls=c("ADN103M2","ADN103M3","ADN103M4","ADN103M5","ADN103M6"))
# exclusion des individus ayant un trop fort pi-hat moyen
system(paste("plink --noweb --bfile ",bfile,"--remove exclude_indiv_meanpi-hat.txt --make-bed --out ",bfile,sep=""))
# exclusion d'individu par rapport au PIHAT
system(paste("plink --noweb --bfile ",bfile,"--remove list_indiv_pihat_sup0.1.txt --make-bed --out ",bfile,sep=""))
# Determination des composantes principales
system(paste("plink --noweb --bfile ",bfile,"--read-genome plink.genome --mds-plot 4 --cluster --out strat",sep=""))
# creation d'un fichier de phenotype à partir des covariables
create_strat_phenofile()
# creation de plots montrant l'association des individus aux composantes, coloration par sex
plot.indiv.components(nb_comp=4, famfile="forStrat2.fam", phenofile="pheno_strat.txt", outputfile="qqplot_components_col_by_sex", color="sex")
# creation de plots montrant l'association des individus aux composantes coloration par status
plot.indiv.components(nb_comp=4, famfile="forStrat2.fam", phenofile="pheno_strat.txt", outputfile="qqplot_components_col_by_status", color="status")
# creation de plots montrant l'association d'individus aux composantes en fonction de l'etude desquels ils sont tires
plot.indiv.components.orig(nb_comp=4, all_famfile="all.fam", casesctrls_famfile="forStrat2.fam", ctrls_famfiles=c("oba.fam","canada.fam"), phenofile="pheno_strat.txt",outputfile="qqplot_components_orig")
# creation d'un unique plot avec une ligne horizontale ou verticale
plot.indiv.components.abline(comp_number=c(1,3),famfile="forStrat2.fam", phenofile="pheno_strat.txt",color="status",Hline=0.02,Vline=-0.04)
# creation de la liste d'individus à eliminer (en fonction des associations aux composantes)
pheno <- read.table(file="pheno_strat.txt",h=T)
c1 <- unique(pheno[pheno$C1 > 0.02,c("FID","IID")])
c2 <- unique(pheno[pheno$C2 < -0.04,c("FID","IID")])
write.table(c1,file="indiv_c1.txt",quote=F,row.names=F) 
write.table(c2,file="indiv_c2.txt",quote=F,row.names=F)
out <- "forStrat3"
system(paste("plink --noweb --bfile ",bfile," --remove indiv_c1.txt --make-bed --out ",out,sep=""))
bfile <- out
system(paste("plink --noweb --bfile ",bfile," --remove indiv_c2.txt --make-bed --out ",bfile,sep=""))
# Determination des composantes principales
system(paste("plink --noweb --bfile ",bfile,"--read-genome plink.genome --mds-plot 4 --cluster --out strat2"sep=""))
# creation d'un fichier de phenotype à partir des covariables 
create_strat_phenofile(mdsfile="strat2.mds",outputfile="pheno_strat2.txt")
# creation des plots
plot.indiv.components(nb_comp=4, phenofile="pheno_strat2.txt", famfile="forStrat3.fam", outputfile="qqplot_components_col_by_sex_cor", color="sex")
plot.indiv.components(nb_comp=4, phenofile="pheno_strat2.txt", famfile="forStrat3.fam", outputfile="qqplot_components_col_by_status_cor", color="status")
# analyse d'association avec ces covariables 
system(paste("plink --noweb --bfile ",bfile," --pheno pheno_strat.txt --all-pheno --assoc --out pheno",sep=""))
# creation d'un tableau recapitulatif des statistiques par rapport aux composantes
summary.qassoc(nb_comp=4, assocfile="pheno", outputfile="summary-qassoc.txt")
# visualisation graphique des resultats
qqplot.components(nb_comp=4, IC=0.95, summaryfile="summary-qassoc.txt", outputfile_sq="qqplot_statsq_assoc_components", outputfile_p="qqplot_pval_assoc_components")


## une analyse de plaque
# creation des fichiers cluster - a faire pour les cas et les controles 
create_cluster_file(inbfile="../final2", plate_file="Plaks_N.csv", num_Pdel=11, status=1)
create_cluster_file(inbfile="../final2", plate_file="Plaks_N.csv", num_Pdel=11, status=2)
# creation du fichier de phenotype - a faire pour les cas et les controles 
create_pheno_file(cluster_file="cluster-ctrls.dat", status=1) 
create_pheno_file(cluster_file="cluster-cases.dat", status=2)

# creation des fichiers d'association sous plink
association()

#creation d'un tableau d'informations sur les chi² et les p-values
info_chisq_pval(p.names=c("P2","P3","P4","P5","P6","P7","P8","P9"), status=1)
info_chisq_pval(p.names=c("P21","P22","P23","P24","P25"), status=2)

# creation d'un QQplot par plaque - a faire pour les cas et pour les controles 
qqplot.p(p.names=c("P2","P3","P4","P5","P6","P7","P8","P9"), nb_by_graph=2, IC=0.95, status=1, outputfile="qqplot_plaque_ctrls")
qqplot.p(p.names=c("P21","P22","P23","P24","P25"), nb_by_graph=2, IC=0.95, status=2, outputfile="qqplot_plaque_cases")

# suppresion des SNPs < 5e-08 pour les contrôles
ct <- read.table(file="p-values-controls-first.txt",h=T)
snp <- ct[ct$P2 < 5e-8,"SNP"]
snp2 <- ct[ct$P4 < 5e-8,"SNP"]
snp3 <- ct[ct$P7 < 5e-8,"SNP"]
snps <- c(as.character(snp),as.character(snp2),as.character(snp3))
write.table(snps,file="del_snps_plaques.txt",quote=F,row.names=F)

qqplot.p(p.names=c("P2","P3","P4","P5","P6","P7","P8","P9"), nb_by_graph=2, IC=0.95, status=1, outputfile="qqplot_plaque_ctrls_clear")


## une analyse finale
# elimination des SNPs en fonction des F_U et F_A 
del_freq(inbfile="../pvmcont_new", th_freq=0.05, outputfile="SNPoutfreq0.05.txt")
system(paste("plink --noweb --bfile ../pvmcont_new --exclude SNPoutfreq0.05.txt --make-bed --out pvm0.05",sep=""))

# recuperation de l'ensemble des donnees 
create_summary_data(inbfile="pvm0.05", del_AT_GC=F, del_heterosome=F, outputfile="summary_data.txt")

#creation d'une liste de SNPs à eliminer suivant Hardy-Weinberg, pourcentage de missing et maf 
snp_selection(HW=0.0001, miss=0.05, maf=0.1,summary_file="summary_data.txt")
system(paste("plink --noweb --bfile pvm0.05 --exclude exclusionSNP.txt --make-bed --out pvm0.05select",sep=""))

# elimination des individus ayant un trop fort mean_pihat (provient de la stratification)
system(paste("plink --noweb --bfile pvm0.05select --remove exclude_indiv_meanpi-hat.txt --make-bed --out final",sep=""))

# elimination des individus par rapport à leur pihat (provient de la stratification)
system(paste("plink --noweb --bfile final --remove list-indiv-pihat-sup0.1.txt --make-bed --out final2",sep=""))

# elimination des individus par rapport aux composantes (C1 : > 0.028 et C2 : < -0.04) (provient de la stratification)
write.table(indiv,file="indiv_C1.txt",quote=F,row.names=F)#dans communs/out
system(paste("plink --noweb --bfile final2 --remove indiv_C1.txt --make-bed --out final2",sep=""))
system(paste("plink --noweb --bfile final2 --remove indiv_c2.txt --make-bed --out final2",sep=""))

# elimination des SNPs par rapport aux plaques (provient de l'analyse par plaque)
system(paste("plink --noweb --bfile final2 --exclude analyse_plaques/del_snps_plaques.txt --make-bed --out final3",sep=""))


######################## ANALYSE D'ASSOCIATION ############################

# association regression logistique finale
system(paste("plink --noweb --bfile final3 --logistic --covar pheno_strat10.txt",sep=""))

# creation des plots de résultats de l'analyse d'association
plots.association.analysis(associationfile="plink.assoc.logistic", qqplot_title="QQPlot", mhtplot_title="Manhattan plot", out_qqplotfile="qqplot", out_mhtplotfile="mhtplot")

